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Abstract. We consider bound geodesic orbits of test masses in the exterior 
gravitational field of a rotating astronomical source whose proper angular momentum 
varies linearly with time. The linear perturbation approach of Lense and Thirring 
is herein extended to the nonstationary case. In particular, we investigate the 
^ ' instability of Lense-Thirring precessing orbits due to the slow temporal variation of 

the gravitomagnetic field of the source. 
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Nine decades ago, Lense and Thirring considered the motion of a free test mass in 
the stationary exterior gravitational field of a rotating astronomical source within the 
OO , framework of general relativity [H |2] . They treated the influence of the gravitomagnetic 
field on the particle orbit via the Lagrange planetary equations and showed by means 



of linear perturbation theory that — when averaged over the fast Keplerian motion — the 
orbit keeps its shape and slowly precesses. Specifically, they found that the semimajor 
axis a of the osculating ellipse does not change, while its eccentricity e and orbital 
inclination i contain periodic terms that vanish on the average; moreover, the osculating 
ellipse precesses. This occurs both within the orbital plane of the osculating ellipse and 
without, as the orbital plane precesses about the axis of rotation of the central source. 
The frequency of both precessions can be described by 

2G 

"^^^ = Q3(i_e2)3/2 [Jo - 3(Jo ■ (1) 

Here, Jq is the constant angular momentum of the source, n is a unit vector parallel 
to the orbital angular momentum of the osculating ellipse and Jq • n = cosi. Thus the 
Runge-Lenz vector and the orbital angular momentum vector of the osculating ellipse 
both precess with the Lense-Thirring frequency ([I]). 

Astronomical bodies in general rotate; however, the magnitude of the proper 
angular momentum is seldom constant. In two recent papers [H H], the gravitational 
physics around a rotating central source whose spin angular momentum vector is fixed 
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in space but varies linearly in time has been explored. In particular, it has been shown 
in [3] that sufficiently far from such a source, the spacetime metric is given by 

ds"^ = -c^ ^1 - 2^^ dt^ - ^(A ■ d^)dt +{^ + Sijdx'dx\ (2) 

where 

^ GM ^ Gjmxx 

$ = , A = 3 

r c r'^ 

are the gravitoelectric and gravitomagnetic potentials, respectively. Here, r = |x|, M is 
the mass of the source and its angular momentum is given by 

3{t) = (Jo + Jit)z; (4) 

moreover, $ < and |A| < c^. Thus r > GM/c^, r > J/{Mc) and all 0(c"^) 
contributions to the metric tensor have been neglected in this linear post-Newtonian 
approach to general relativity. As explained in [3l H] , we simply ignore the processes by 
which the variation of angular momentum is turned on and off and assume that equation 
(jlj) holds throughout the temporal interval of interest; furthermore, all radiative effects 
are neglected. 

The motion of a free test particle in the gravitational field of the source is given by 

the geodesic equation in the spacetime with metric ([2]). This equation, as shown in [1], 

can be written in its reduced form 

dv GM^ GM , , , . , 2G j X X 2 

— + = ^ 4(x ■ v)v - v^^] + —— V X B 

dt r'^ c r'^ c r'^ c 

-^[J-(-xv)](x-v)v, (5) 

where an overdot represents differentiation with respect to time t and B = V x A is 
the gravitomagnetic field given by 

G{Jo + Jit) 2m r(i\ 

B = = (32X — rz). fo) 

cr° 

The right-hand side of (JSj) contains all of the linear post-Newtonian contributions 
that arise from potentials given in ([3]). It turns out, however, that in a general treatment 
to 0(c~^), the nonlinear gravitoelectric term 4G^M^x/(c^r^), which is quadratic in $ 
and hence absent in our linear treatment, should also be taken into account. In the 
present work, we explore further the influence of the temporal variation of J on motion 
around a central rotating source to first post-Newtonian order, namely, 0(c~^). Thus 
instead of ([5]), we consider 

dv GMx _ 

M^ — -^' « 

„ GM, , , , , AG^M"^ 2G j X X 2 

[4(x ■ v)v - v^^] + X + — — vxB. (8) 



c^r^ c^r^ c 

As demonstrated in [3l H], equation ([2]) represents the metric of a nonstationary 
linearized Kerr spacetime. The geodesic equation in Kerr spacetime is completely 
integrable [5|; more recent results are contained, for instance, in [6j and references 
therein. 
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2. Instability of spherical orbits 

To gain insight into the nature of allowed orbits, it proves useful to study first certain 
simple configurations. For instance, it has been shown in [1] that circular equatorial 
orbits are unstable due to the temporal variation of J. We therefore look for spherical 
orbits in this section. It will turn out that the post-Newtonian equations of motion ([71) 
and ([8]) do not allow spherical orbits unless J = 0. For J = 0, the spherical orbits 
can be concisely characterized as describing a post-Keplerian circular orbit undergoing 
Lense-Thirring precession. We will then investigate the instability of spherical orbits 
under the slow temporal variation of J. 



2.1. Equations of motion in spherical coordinates 

To investigate spherical orbits, it is convenient to write ([7]) and ([8]) in spherical polar 
coordinates (r, 9, 0) such that x = r sin 9 cos (p, y = r sin 9 sin (p and z = r cos 9. In terms 
of the corresponding unit vectors f , and (f), one can write x = rf , 

X = V = rf + r90 + r0 sin 9 <p, (9) 
X =(f — r9^ — sin^ 9)f + {r9 + 2r9 — sin 9 cos 9)6 
+ (r0 sin 9 + 2r0 sin 9 + 2r9(j) cos 9) (p. 
Moreover, it follows from ([8]) that F = F^f + FgO + F^<p, where 
^ GM f ^ 2 ^GM\ 2GJ(t) • 2 . 

Fr = ^Ur^-v^ + ^—]+ , y ^sin'g, (11] 



•■ A2 12 ■ 2n GM GMf^., 2 4G'M\ 2GJ{t)^.,,^ 
r-r9^ - rcj)^ sm^ 9 + = — - Ar^ - + + — ^4^0sm^ 9, (14) 



dt 



21 ■ 2n '^GJ{t) . 2^ 

r sm 9 sm 9 

c^r 



C2 



(10) 



^ 4GM ■ 4GJ(t) ■ ^ ^ 

Fe = -^r9 ^^sin^cos^, 12 

= ^5 — r0sm6' + — — sm6' ^-^(rsm6' - 2r9cos9). (13) 

Thus the post-Newtonian equations of motion are 



r^' + 2f^ — sin 9 cos ^ = — - — r9 sin 9 cos 6', (15) 

AGM • 2Gj 
rcj) sin 6' + 2r0 sin ^ + 2r^0 cos ^ = — - — r0 sin ^ H — — — - sin ^ 

~'^^^{rsm9-2r9cos9). (16) 
It is interesting to note that (flGil may be written as 



rr0 sm 9. (17) 



Thus an integral of the motion is obtained if the right-hand side of (fTTjl vanishes. We 
therefore look for spherical orbits (r, 9, 0) = (p, "i?, v?) such that p is a constant. 
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2.2. Spherical orbits 

It follows from ( fTTl) that for a spherical orbit {p,{},ip), 
C 2GJ{t) 

¥'=-^ + ^T^' (18) 

where C is a constant of integration. Substituting (fT8|) in (fHj) and (ITSl) . and keeping 
terms only up to 0(c~^), we find 

sm w P \ c^P / c^P 
''-^^-O. (20) 



and 

sin^ ^9 

respectively. These equations are compatible only if J = 0. Thus let J = Jq and define 
the positive constant f2 > such that 

tf.^fi-^^^-^^a (21) 



p-^ \ c^p J c^p^ 
To find the motion in we note that f|T9|) can be written as 



={Sl' -C')-n'cos'^, (22) 

which has solutions only for > C^; in this case, we have 

cos^9 = asin(fit + /?), (23) 
where a is given by 

C2 = fi2(i_«2) (24) 

and /3 is a constant phase. It is possible to let C = Qcosi and a = ±sini, where z is a 
constant ( "inclination" ) angle. Then the solution of ( ITSj) may be expressed as 

2GJ t 

ip{t) = ° + tan"^ [cos i tan(fit + p)] + (fo, (25) 
where ipo is an integration constant. Let us note that with 

where ujk is the Keplerian frequency, we have 

SGJpCosi 

^ = ^ ^1 — • 27 

Only positive square roots are considered throughout. The explicit solution for spherical 
orbits in the first post-Newtonian approximation is contained in ( 123|) and ( l25l) . Spherical 
orbits in Kerr spacetime have been considered in [7]. 

It is interesting to note that when Jq = in ( |23ll and ( |25l) . a spherical orbit 
simply reduces to a circular orbit of radius p and frequency uj about mass M such that 
the orbital angular momentum vector makes an angle of i with the z axis. Thus one 
can characterize the spherical orbits under consideration as circular orbits in the post- 
Newtonian gravitational field of mass M that undergo Lense-Thirring precession due to 
the presence of the constant angular momentum of the source Jq. 
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2.3. Perturbed spherical orbits 

We now let J = Jq + Jit with Ji 7^ and consider the solution of equations ( fT4l) -( fT7l) 
to linear order of perturbation beyond an arbitrary spherical orbit. Thus let 

r = p(l + /), e = ^ + g, (j) = ip + h, (28) 

where fit), g{t) and h{t) are all of 0(c~^) and depend linearly upon Ji. Substituting 
( I28l) in the equations of motion, we find 

f - 2dg - 2Ch - Suy - 2^g = (29) 

AC J Ct 

g-2^ coi{2^)g + 2^f - 2C cot = -i— cot t?, (30) 

c^P'^ 



^ + -^9 + -2-J = -^ + —i- (31) 



2 •■ 2C , 2^71^ CRo 
C sm V c^p'^ sm w 

Here, C = ucosi and Rq is a dimensionless integration constant. While equations ( l29l) 
and (l30l) have been respectively obtained directly from (fl^ and (ITS!) . equation (!3T!) 
is obtained from ( JT71) . since the right-hand side of this equation vanishes under our 
perturbation conditions. Substituting for h in fl29l) and (l30l) using fl3T|) leads to 

■r f 2 \ r ra a . QG JlCt 2RqC^ , , 

\ sm V j c^p'^ sm it 

g + 2^2 cot ^ - coi{2^)]g + 4^/ + 2^f = 2Rq^. (33) 
Inspection of equations ( l32l) and ( l33i) reveals that they have the following solutions: 

fit) = ^^-^Kt + 2i?o, (34) 

^7(t) = - {^-^j^^Kt' + 3i?ot - 60) ^, (35) 

where Go is an integration constant. With these formulas for / and g, ( l3Ti) can be 
simply integrated and the result is 

, , , GJit^ C /6Jicosz 9 ^ \ 

= -h ^ 1,2 + 3^ot - eo + $0, 36 

where $0 is another constant of integration. 

It is possible to assume that i?o = ©o = "^"0 = without any loss in generality. To 
see this, let Ji = 0; then, it is straightforward to demonstrate, using equations fl28|) and 
flM|) - fl5B]) . that the resulting orbit (r, 6*, 0) is simply a new spherical orbit of constant 
radius p(l + 2_Ro)- K follows that one can set the dimensionless constants Rq, Gq and $0 
equal to zero with no loss in generality. Then, the perturbation (/, g, h) due to Ji 7^ 
is clearly secular leading to the instability of spherical orbits; in fact, these orbits spiral 
outward for Ji cosi > and inward for Ji cosi < 0. These results provide independent 
confirmation of the some of the conclusions regarding perturbations of Keplerian orbits 
presented in since a spherical orbit in the present treatment is simply a Lense- 
Thirring precessing circular orbit. A detailed investigation verifies this correspondence 
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for the secular terms. We note, in this connection, that for Ji = 0, the perturbed 
circular Keplerian orbit in [3] differs from a spherical orbit by a harmonic term in the 
radial perturbation — see equation (71) of [3]. 

Finally, let us remark that the speed of the motion along the perturbed orbit is 
given by 

2 2 SGJiCosi 
V =Vf;, 5 ujRt, (37) 

where Vq is the speed of the unperturbed spherical orbit 

^o = P^(l-^^cosz^ . (38) 

It follows from fl37|) that the motion is slower for Ji cosi > and faster for Ji cosz < 0. 
This confirms one of the main results of [4j by a completely different analysis. On these 
physical grounds, we expect that our special solution contains the dominant secular 
terms of the general solution of the linear perturbation equations. 



3. Extension of the Lense-Thirring approach 

We turn now to the gravitomagnetic perturbations of Keplerian ellipses, first treated in 
general by Lense and Thirring [H [2] • In this section, we extend their analysis to take 
into account the temporal variation of the gravitomagnetic field. 

The method that we employ here is not a generalization of the approach developed 
in the previous section for circular orbits. That is, rather than directly perturbing a 
Lense-Thirring precessing Keplerian ellipse, we follow, for the sake of simplicity, the 
linear perturbation method of Lense and Thirring [H [2] in the case of J(t) = Jq + Jit. 

It is possible to express the three second-order equations of motion in terms 
of six first-order Lagrange planetary equations p]. If at any instant of time t the 
perturbing force F is turned off, the test particle follows an osculating Keplerian ellipse 
about the central source. Thus instead of the position and velocity of the particle at 
time t, the state of the particle can be equally well characterized by the six orbital 
elements of the osculating ellipse at time t. The motion can therefore be described 
in terms of the evolution of the parameters of the instantaneous osculating ellipse. 
These parameters can be chosen in various ways; we employ the Delaunay action-angle 
variables {L,G, H, £,g,h) given by [9| 

L = G = [GMa{l - e^)]^/\ H = Gcosi, (39) 

i = u — e sin u, q = argument of the pericenter, 

^ ^ ^ (40) 

h = longitude of the ascending node. 

Here, a is the semimajor axis of the osculating ellipse, e is its eccentricity, i is the orbital 
inclination, u is the eccentric anomaly and £ is the mean anomaly. The quantity G is 
the magnitude of the specific orbital angular momentum vector G = x x v, while its 



a(l — e cosm), (41) 
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2;-component is denoted by H. Moreover, the radial position of the test particle along 
the osculating ellipse is given by 

a(l-e2) 

r = r, 

1 + e cos V 

where v is the true anomaly. 

The equations of motion are 

dL 

— = ^iFrCsinv + Fs(l + ecosv)], (42) 
dt G 

dG ^ 

* = (43) 

— — = r[F5 cosi — F„ sini cos('i; + ^)], (44) 
d£ r 

— = ojk H ^[-F'r(— 2e + COS'S + ecos^f) — Fs{2 + ecos-u) sint>], (45) 

dt e 

'^^^ . • • , -A , i^-e'Y^' ( ^ ~^p2 + ecosS. \ 

— = — cot I sm( t> + g) ^ —rr cos v + r « r sm f , ( 4o ) 

dt G uxcie \ 1 + e cos v J 

dh ^ rF„ sin(ti + g) 

dt ~ G sini ' ^ ^ 

where the Keplerian frequency is given by ujk = (GMY^"^ / . Here, the perturbing 
force is given in terms of its radial, sideways and normal components, 

F = F,r + F,s + F„n. (48) 

That is, r is the radial unit vector as before, s = fi x f and fi is given by G = Gfi; 
hence, r and s are in the plane of the osculating ellipse, while fi is normal to it. 

To find the components of the perturbing force, it is useful to recall that the position 
vector of the test particle along the osculating ellipse has components 

X = r [cos h cos{v + g) — sin h cos isin{v + g)], (49) 
y = r [sin h cos{v + g) + cos h cos ism{v + g)], (50) 
2; = rsinzsin(?) + ^), (51) 

while its velocity is given by 

v = rr+— s, r = GM — (52) 
r G 

Furthermore, with respect to the background (x, y, z) coordinate system 

fi = (sin /i sin i, — cos /i sin i, cosi), (53) 
s = (— cos h sin({) + a) — sin h cosi cos({; + q), 

— sin /isin(t> + g) + cos h cosz cos(f + g), sini cos(f + g)). 

The perturbing force (IHl) contains small relativistic terms and in our perturbation 
scheme the influence of these perturbing accelerations are simply additive. Therefore, 
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we will ignore the post-Newtonian gravitoelectric terms in F and concentrate instead 
on the gravitomagnetic terms, namely, 

^, 2Gjxx 2 

F' = -^—^ V X B. 55 

We find that 

^, 2G( Jo + J,t)H 
= — , 

GMe sin v 



. (56) 



, 2G cos i 
r ^ = 



^2 2 



Ji - (Jo + Jit)- 



(57) 



rG 

, 2Gsini f _^ . GM ~ , _x /~ 

-TV = ^r^\ Ji cosff + g) + (Jo + Ji^) — ^ esin^f - (2 + Secost;) sin(f + ^f) >. (58) 

With these perturbing functions, we can simply integrate (H2l) - fj47j) following the linear 
perturbation approach adopted by Lense and Thirring [U [2]. That is, we regard the 
orbital elements appearing on the right-hand side of the equations of motion as constants 
and employ dv/dt = G/r"^ for the osculating ellipse at time t; in fact, v{t) and t{v) are 
obtained from 

dv U)k / ~\2 

^=(13^(1 + ^^°^^) • (59) 

Assuming, for the sake of simplicity, that t = at 5 = 0, we find 

5 

V = uJxt + 2esinci;i^t + -e^ sin2ct;i^t + 0(e'^), (60) 
3 

ujKt = V - 2esmv + -e^ sin2t; + O(e^). (61) 

The integration of the equations of motion is now straightforward; we recover the 
Lense-Thirring results for Jq and find new terms proportional to Ji. The latter 
terms are generally secular, of course, and render the Lense-Thirring precessions time- 
dependent. In connection with orbital instability, the rest of this section is devoted to 
the secular variation of (a,e, z), which remained on average unchanged in the Lense- 
Thirring treatment. 

It follows from <^ that 

da 4GJia(l — e^)^/^ cosi 

17 = 2 3 ' (62) 

which can be integrated using (159|) and the result is 

4Jiacos2 . . 

= Mc2(l-e2) + esm^;). (63) 

Thus the semimajor axis and hence the Newtonian energy of the osculating ellipse will 
have a secular variation in time depending upon the sign of Ji cosi, so that essentially 
all orbits are unstable for Ji ^ 0. Similarly, from = GMa{l - e^), (jlJ]) and (jlS]), we 
find that the change in eccentricity is given by 

2ci;i^COS2 r X A / ~N T A 2JiCOS2 . ^, 

= Mc^(l-e^)V2 [-'^°^(^"^^) + ^1^^] + -J^^^^ (6^) 
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where 

I=ftWA.MX. f2co.A + e(l+co.^A)^^ 

Jo Jo 1 + e cos A 

These integrals contain secular terms. For instance, relations (!60!) and (!6T]) may be used 
to show that for 5 = at t = 0, 

ukI = sint-tcost + |(sin2t-2tcos2t) + 0(e^), (66) 

where t = uxt. Integral /' can be evaluated exactly (see, for instance, formula 2.553 3 
on page 148 of [ID]); it can also be expressed in powers of e as 

/' = 2sin{} + ^(2{}-sin25) + 0(e2). (67) 

For an orbit with initial eccentricity much less than unity, the secular term that is 
independent of the eccentricity in Ae is given by the variation of 
^Jicosi 

-2 uxtcosujKt. (68) 

Let us also note that Ae and Aa are both proportional to cosi, so that an osculating 
polar orbit tends to preserve its shape. 

Finally, equations (H3l) and (jHj) can be used to show that 

f^ = 'Scos{v + g), (69) 

where is given by (!58l) . A simple integration reveals that 

Ji sin i , ^ uk sini ^ . , . 

where 

Jq 1 + e cos A 

n [ 3 1 

J' = (Jo + Jit) (1 + -ecosA) sin(2A + 2^) - esin^cos(A + ^) rfA. (72) 
Jo L 2 J 

Here, J" can be evaluated exactly [I^; moreover, ^7' can be expressed in powers of e 
using (l60l) and (161]) . We note that Ai is proportional to sini, so that an equatorial orbit 
stays in the equatorial plane. In equations (jTOl) and (|72l) . the part proportional to Jo 
contains only harmonic terms, while the part proportional to Ji contains secular terms. 
For e <^ 1, the secular term that is independent of the eccentricity in (1701) is given by 
the variation of 

Jisinz 2 _ 

- 2 ujKt cos [ujKt + g). (73) 

To further illustrate orbital instability, we integrate numerically the equations of 
motion for an initially eccentric Keplerian ellipse. In practice, it turns out to be simpler 
to integrate equation ([7]) — where F is now replaced by F' given in fl55l) — with initial 
position and velocity given by ( H9l) -(l52 l) . Moreover, we employ the equations of motion 
in dimensionless form, so that all lengths are given in units of Oq, the initial semimajor 
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Figure 1. Plots of a/ao, e, i and v of the osculating ellipse, respectively from the top 
panel down, versus lo^t for the parameters given in the last paragraph of section 3. 
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axis, while time is given in units of l/cu^, the inverse of the initial Keplerian frequency. 
The equations of motion then depend on two dimensionless parameters 

2Jouj% 2Ji 
^° = -M^^ = M^- ^^^^ 

For the numerical results illustrated in figure 1, we choose 6q = 10^^. This 
corresponds approximately to an initial orbit of semimajor axis oq = 40 km around 
a neutron star of mass M ^ 2Mq and radius ~ 10 km with a proper rotation period of 
a millisecond. Furthermore, we choose 6i = — 10~®, so that after about 1600 Keplerian 
periods, the angular momentum of the source decreases to zero. The integration is 
carried out for u^t : 20000 such that J(t) : Jq — Jq. This relatively rapid 
decrease of angular momentum has been adopted here for the sake of illustration; in 
fact, neutron stars generally lose angular momentum very slowly due to electromagnetic 
braking torques. In figure 1, we plot a/ao, e, i and the speed of the motion v versus 
uj^t. For the initial conditions at t = and v = 0, we choose ho = vr/6, go = 7r/3 and 
io = 7r/4. The initial eccentricity is chosen to be Cq = 0.5; based on our numerical work, 
similar results are expected for other initial eccentricities. The simple linear behavior 
of a{t) depicted in the first panel of figure 1 can be obtained from (|63|) : that is, for 
uj^t ^ 1, equation fl63l) implies that 

/ 2(5icosio A /vp-N 
"0 ( 1 + — —^xt j , (75) 



X - Co 

in agreement with our numerical results. The second and third panels of figure 1 depict 
the oscillatory character of the eccentricity and the inclination angle, respectively, as 
the angular momentum of the source monotonically decreases from Jq to — Jq. The 
midpoint of integration when J = is a prominent feature of these graphs. In fact, the 
amplitudes of the quasi-periodic oscillations appear to be proportional to J. Moreover, 
the inclination angle i tends to oscillate toward the angular momentum vector of the 
source J. As indicated by our extensive numerical work, and is evident from the last 
three panels of figure 1, we have to expect complexity in the details of the motion, which 
is probably chaotic; therefore, only the overall trends are meaningful here. The last 
panel confirms the expectation that for Ji cosi < 0, f on the average has an increasing 
trend with time, in general agreement with equation ( l37|) for the circular case. That 
is, for Jicosi > (Jicosi < 0), the orbit generally tends to spiral outward (inward) 
accompanied by a corresponding decrease (increase) in its average speed. 



4. Discussion 



We have studied the instability of bound Keplerian orbits induced by a time-varying 
gravitomagnetic field in the post-Newtonian approximation. Circular and elliptical 
orbits have been treated separately in sections [2] and [3l respectively. The results are 
expected to be of interest in the study of variable collapsed astrophysical systems. 
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